% Main_Function based on GPU Computing
% It is noted that matrix operation need to be performed in the divided blocks to avoid video card's memory overflow.
% 'gpu**' shows gpu data;
% W.Q. 2020.09.02

% w: the external excitation frequency
% Block_Size: the size of the Block to avoid video card's memory overflow. 
 
function [X0,X1,X2,Y0,Y1,Y2] = Grid_Method(w,X0Hat,X1Hat,X2Hat,Y0Hat,Y1Hat,Y2Hat,Block_Size,tolerance,Data_File_Number)

	g = gpuDevice;
	gpuThreshold_Value=gpuArray(tolerance);                                               % the threshold value of period condition for nonlinear force
	% gpuThreshold_Ratio1=gpuArray(0.1);                                              % the threshold ratio 1 of period condition for nonlinear force
	% gpuThreshold_Ratio2=gpuArray(1.1);                                              % the threshold ratio 2 of period condition for nonlinear force
	X = [ X0Hat;X1Hat;X2Hat;];
	Y = [ Y0Hat;Y1Hat;Y2Hat;];

% L, DL, PINV L
	NP = 1024;ND = 1;NH = 1;v = 1;
	[ts,L,DL,PL] = FFT_Operator(NP,ND,NH,v,w);  
	gpuL=gpuArray(L);                                                               % L of CPU -> L of GPU
	gpuDL=gpuArray(DL);                                                             % L of CPU -> L of GPU
%	gpuPL=gpuArray(PL);                                                             % PL of CPU -> PL of GPU
% 	gpuX=gpuArray(X);                                                               % X of CPU -> X of GPU
% 	gpuY=gpuArray(Y);                                                               % Y of CPU -> Y of GPU
	gpuw=gpuArray(w);                                                               % w of CPU -> w of GPU
	gputs=gpuArray(ts);                                                             % ts of CPU -> ts of GPU
	gpuIndicator_1=zeros(1,size(X,2),'gpuArray');                                % No.1 indicator of GPU
	gpuIndicator_2=zeros(1,size(Y,2),'gpuArray');                                % No.2 indicator of GPU
% gpuIndicator_3=zeros(1,size(X,2),'gpuArray');                                % No.3 indicator of GPU
% gpuIndicator_4=zeros(1,size(Y,2),'gpuArray');                                % No.4 indicator of GPU
	index_Block = buffer(1:size(X,2),Block_Size);
	for iBlock = 1 : size(index_Block,2)
		gpuX=X(:,index_Block(index_Block(:,iBlock)>0,iBlock));
		gpuY=Y(:,index_Block(index_Block(:,iBlock)>0,iBlock));
		gpuxs=gpuL*gpuX;             % This command is equivalent to 'index_Block(find(index_Block(:,i_Block)>0),iBlock)'
		gpuvxs=gpuDL*gpuX;
		gpuys=gpuL*gpuY;
		gpuvys=gpuDL*gpuY;

		[ gpuFxs,gpuFys ] = Bearingforce(gputs,gpuw,gpuxs,gpuys,gpuvxs,gpuvys);
		gpuDet=gputs(2)-gputs(1);
		gpuIndicator_1(1,index_Block(index_Block(:,iBlock)>0,iBlock)) = sum(gpuDet.*gpuFxs(2:end,:))-mean(gpuFxs).*gputs(end);          % Matrix version
		gpuIndicator_2(1,index_Block(index_Block(:,iBlock)>0,iBlock)) = sum(gpuDet.*gpuFys(2:end,:))-mean(gpuFys).*gputs(end);          % Matrix version
		% gpuIndicator_3(1,index_Block(index_Block(:,iBlock)>0,iBlock)) = abs(gpuFxs(1,:)-gpuFxs(end,:))./abs(gpuFxs(end,:));          % Matrix version
		% gpuIndicator_4(1,index_Block(index_Block(:,iBlock)>0,iBlock)) = abs(gpuFys(1,:)-gpuFys(end,:))./abs(gpuFys(end,:));          % Matrix version    
	end
	gpuIndex_SET1=find(abs(gpuIndicator_1)<gpuThreshold_Value);
	gpuIndex_SET2=find(abs(gpuIndicator_2)<gpuThreshold_Value);
	% gpuIndex_SET3=find(abs(gpuIndicator_3)<gpuThreshold_Ratio1 | (abs(gpuIndicator_3)<gpuThreshold_Ratio2 & abs(gpuIndicator_3)>1) );
	% gpuIndex_SET4=find(abs(gpuIndicator_4)<gpuThreshold_Ratio1 | (abs(gpuIndicator_4)<gpuThreshold_Ratio2 & abs(gpuIndicator_4)>1) );
	gpuIndex_SET=intersect(gpuIndex_SET1,gpuIndex_SET2);
	Index_SET = gather(gpuIndex_SET);                                             % Index_SET of GPU -> Index_SET of CPU
	Reduced_X = X(:,Index_SET);
	Reduced_Y = Y(:,Index_SET);
	reset(g);
  % Debug
	X0=Reduced_X(1,:)';
	X1=Reduced_X(2,:)';
	X2=Reduced_X(3,:)';
	Y0=Reduced_Y(1,:)';
	Y1=Reduced_Y(2,:)';
	Y2=Reduced_Y(3,:)';  
	save(['Samples',num2str(Data_File_Number),'.mat'],'X0','X1','X2','Y0','Y1','Y2');
end